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ABSTRACT 

We have developed a modified form of the equations of smoothed particle magnetohy¬ 
drodynamics which are stable in the presence of very steep density gradients. Using 
this formalism, we have performed simulations of the collapse of magnetised molecu¬ 
lar cloud cores to form protostars and drive outflows. Our stable formalism allows for 
smaller sink particles (< 5 AU) than used previously and the investigation of the effect 
of varying the angle, '^9, between the initial field axis and the rotation axis. The nature 
of the outflows depends strongly on this angle: jet-like outflows are not produced at all 
when if > 30°, and a collimated outflow is not sustained when if > 10°. No substantial 
outflows of any kind are produced when if > 60°. This may place constraints on the 
geometry of the magnetic field in molecular clouds where bipolar outflows are seen. 

Key words: accretion, accretion discs - MHD - stars: formation - stars: winds, 
outflows. 


1 INTRODUCTION 

Magnetic fields are one of the most important forces influ¬ 
encing the formation of protostars and may resolve several 
questions about the formation of stars that are left unan¬ 
swered by purely hydrodynamic theories. That the molecu¬ 
lar clouds that ultimately produce protostars are magnetised 
is well known (Crutcher et al. 1993; Crutcher 2012), and this 
can be confirmed with observations of Herbig-Haro objects 
with distinctive bipolar outflows, which must have a mag¬ 
netic origin. Additionally, these outflows may help explain 
the difference between the angular momentum observed in 
molecular cloud cores and that of the resultant stars. This 
is, however, dependent on the ability of the protostar to pro¬ 
duce a strong outflow - if at certain angles this is impossible 
this may place constraints on the initial field geometry. Re¬ 
cent advances in observational technology have shown that 
this magnetic field structure can be quite complex (Stephens 
et al. 2014) and that the hitherto common assumption that 
field, outflow, and rotation axis are all aligned may be in¬ 
correct (Hull et al. 2013). On very small scales Donati et al. 
(2010) have observed a 20° misalignment between the rota¬ 
tion and field axes of AA Tan. Previous work, e.g. Ciardi 
Sz Hennebelle (2010), using AMR codes has shown that the 
nature and extent of the protostellar outflow is strongly de- 
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pendent on the angle (which we denote with if) between the 
field and rotation axis. 

Smoothed particle hydrodynamics (SPH) methods have 
been applied to many problems related to the formation of 
stars, beginning with the original work of Lucy (1977). These 
hydrodynamic methods have been extended to magnetohy¬ 
drodynamics (MHD), originally by Gingold & Monaghan 
(1977) and Phillips & Monaghan (1985), with limited suc¬ 
cess. Hitherto, such work has been limited by various nu¬ 
merical instabilities, ranging from unphysical pairing (the 
Tensile instability’) of SPH particles (Swegle et al. 1995; 
B0rve et al. 2001) to the production of equally unphysical 
non-solenoidal fields (Tricco & Price 2012). The most recent 
instability, and the one that this paper tackles is that a for¬ 
malism of smoothed particle hydrodynamics incorporating 
fixes to all of the above deficiencies (Price 2012) is unstable 
when ‘small’ sink particles (Bate et al. 1995) are used. 

Most recently. Price et al. (2012) examined the collapse 
of a magnetised molecular cloud core all the way to the for¬ 
mation of the first hydrostatic core (Larson 1969) and Bate 
et al. (2014) have continued to the stellar core. To model the 
evolution significantly beyond protostar formation, sink par¬ 
ticles are a necessary evil since modelling both the magne¬ 
tohydrodynamics of the protostar and also the surrounding 
cloud is computationally unfeasible due to the widely dif¬ 
ferent length and time scales involved. Consequently, some 
way of stabilising the equations of smoothed particle mag- 
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netohydrodynamics (SPMHD) in these cases is essential to 
make progress. These instabilities seem to be only magni¬ 
fied by mis-aligned fields. Previous SPMHD modelling of 
collapsing cores use somewhat large (5 AU) sink particles 
and thus, whilst stable, failed to capture the full range of 
physics. As well as providing a deeper understanding of the 
formation of individual stars, such a formalism could then 
be used in larger, cluster size, simulations similar to Bate 
(2012) but with the addition of magnetic fields. Previous 
cluster-scale simulations performed with SPMHD used the 
Euler potential method which is limited to certain field ge¬ 
ometries (Price & Bate 2008, 2009). Given the huge range 
of density present in following a magnetised cloud collapse, 
a Lagrangian method such as SPH is ideal. Using our modi- 
hed SPMHD method we are able to follow the collapse much 
further than previously, with arbitrarily small sink particles, 
and at a much higher resolution. 

In section 2 we describe our SPMHD formalism, the 
cause of the instability seen in previous work, and the modi¬ 
fications made to eliminate this. Section 3 details our initial 
conditions. We then perform a low resolution test using a 
differentially rotating ‘accretion disc’ and also a collapsing 
magnetised cloud core to demonstrate this modification in 
section 4. Finally, in section 5 we apply this new formalism 
to the collapse of magnetised cloud cores with several dif¬ 
ferent values of ^ and we discuss the effects of varying this 
parameter. 


2 METHOD 

2.1 Standard SPMHD 

As in Price et al. ( 2012 ), we evolve the equations of ideal 
magnetohydrodynamics with the addition of gravity, viz. 

( 1 ) 


- V> 

dt p 

( 2 ) 

= (bW^) v' - B' (v^w^) , 

( 3 ) 

= 47rGp , 

( 4 ) 

with the MHD stress tensor is given by 


+ —( B*B^ - b^-^Bb 

Mo V 2 J 

( 5 ) 

and where, as usual, p, ^ ^ P, (j) represent density, veloc¬ 

ity, the magnetic field strength, the hydrodynamic pressure 
and the gravitational potential, respectively; repeated indi- 
cies imply summation, and 

d d i—i 

dt = m+^^ ■ 

(6) 


We evolve these equations using the method of 
smoothed particle magnetohydrodynamics described in 
Price & Monaghan (2005) and Price (2012) and artificial 
viscosity and resistivity terms based on Riemann solvers 
(Monaghan 1997) with temporal and spatially dependent 
switches. We use the Morris & Monaghan (1997) switch 


for artificial viscosity, with aAv ^ [ 0 . 1 , 1 . 0 ] and the newer 
Tricco & Price (2013) switch for artificial resistivity with 
as G [0.0,1.0]. This differs slightly from Price et al. (2012) 
which used an older resistive switch and constrained an 
to [0.0,0.!]. We soften the gravitational potential using the 
same SPH smoothing kernel as used in the rest of the sim¬ 
ulation (Price & Monaghan 2007). 

All magnetic fields are solenoidal and consequently we 
must maintain a divergence free field. This is not natu¬ 
rally satisfied in SPMHD. Consequently we use the con¬ 
strained hyperbolic divergence cleaning of Tricco & Price 
(2012) (which is based on the earlier Dedner et al. (2002) 
method used in some grid codes). With this, we have an ad¬ 
ditional scalar field, -0 in the induction equation such that: 

^S%lean = -Vy , (7) 

where 

-t-lt (V‘„‘) . (8) 

This removes the unphysical divergence by propagating a 
damped wave through the simulation. Unlike in Bate et al. 
(2014), we do not need to increase the cleaning wave, Cc, 
speed above the magnetosonic speed to maintain stability, 
avoiding a costly decrease in the size of the timesteps. We 
set the damping timescale. 


to be critically damped with a = 0.8, where h is the SPH 
smoothing length. 

As in Price Sz Monaghan (2004), we use a variable 
smoothing length formalism to ensure that computational 
resources are used efficently and that sufficent resolution is 
applied to the complicated areas of the model. The smooth¬ 
ing length and density are solved self-consistently (via the 
Newton-Raphson iterative method) using 



where 77 = 1.2 for the ‘standard’ cubic B-spline kernel (Mon¬ 
aghan 1985), and = 3 is the number of spatial dimensions. 

The simulations were performed using a three- 
dimensional SPH code - originally written by Benz et al. 
(1990) but extensively modified by Bate and his collabo¬ 
rators (Bate 2009) - with the additional modifications de¬ 
tailed in the next section. The code uses a binary tree to 
both hnd neighbours for particles and to calculate gravity. 
A second-order Runge-Kutta-Fehlberg integrator (Fehlberg 
1969) with each particle carrying an individual timestep 
(Bate et al. 1995) was used to evolve the simulation in time. 
Each simulation was run on a single 12-core hyperthreaded 
compute node (i.e. 24 execution threads in total), taking 
between 400 hours of wall time (4,500 core hours) for the 
simpler aligned models to over 550 wall hours (6,600 core 
hours) for the more complicated highly misaligned models. 


2.2 The ‘average h’ method 

As noted earlier, an SPMHD instability exists in regions 
where the density gradient is very large. Since p and the 
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smoothing length, h are related by eq. (10) any density 
gradient will naturally produce an inverse gradient in the 
smoothing length. For conventional SPH, and indeed for 
SPMHD where the gradients are more gentle, this does not 
produce any stability problems. To discover why this is an 
issue, we consider the correction to the tensile instability 
alluded to earlier. 


Separating the MHD stress tensor into two parts (for 
the remainder of this paper we will ignore the gravitational 
potential terms) such that 


gij ^ gij 


isotropic 


+ 5'^'I 


anisotropic 5 


( 11 ) 


where 

( 12 ) 


(13) 

Mo 

we observe that the anisotropic component of the momen¬ 
tum equation (eq. (2)) can be written as 



anis 


= -i— 

P PPO 


(14) 


where we note the important constraint - V • 0 - that 

implies that the term must be zero to ensure a 

solenoidal field. As noted before this is not true in general 
in SPH and an unphysical force along the field lines may 
be produced - this is the Tensile instability’ (Swegle et al. 
1995). Specifically, this results in the stress tensor becom¬ 
ing positive, and hence an attractive force between particles 
being generated. The cleaning described earlier is correcting 
for a different manifestation of the same problem, and hence 
will not help here. The correction proposed by Bprve et al. 
(2001) is to subtract a source term which is exactly equal to 
this unphysical divergence. In our formulation of SPMHD 
we use a symmetric operator for the anisotropic component 
of the momentum equation, viz. 


dt 


'^a I anis,full 


1 ^ 

= -E 

Mo V 


rrih 


+ 


^hpl 


^ViWabiha) 

^ ^aPa 


ViWabihb) 


(15) 


where Wab {h{a,b}) = W (r’ -rl,hia,b}) is the smoothing 
kernel, a and b represent individual SPH particles, N is the 
number of neighbour particles (i.e. particles for which Wab ^ 
0 ), 

O _1 dha^dWabiha) 

“ dpa A dha 

_ ^ hg dWab (hg) 

lypg 4^ dhg 

b 

are terms to take account of gradients in h, is the number 
of spatial dimensions, and all other terms have the usual 
meanings. Consequently, we use a symmetric operator to 
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estimate the magnetic divergence, 

ViBi = —J2mb{ fP^^iWab (ha) 

Mo V \^apl 

^ /-I 

+ ^'^iWab (hb) 

^hpb 

which is then subtracted from the momentum equation such 
that, 

^«ianis = ^^iania.full " (viB^) . (18) 

We use X = 1 as recommended by Tricco & Price (2012). 
Whilst this necessarily makes the equation non-conservative, 
it is much more stable than the value of | recommended by 
Bprve et al. (2001). This means the anisotropic momentum 
equation we evolve is given by 

d 1 ^ 

Comparing this to the SPMHD induction equation, 

is (S)=-nil y- ('■• - “O 

we observe that eq. (20) depends only on hg and eq. (19) only 
upon hb. In situations where the gradients in p are small, this 
does not present any major issues. However, if pa ^ pb then 
ha hb, consequently, it is possible for some particle a to 
have the temporal evolution of its magnetic field evaluated 
over a very small number of neighbours but a force from 
that field interpolated over a large number. This is clearly 
undesirable, and is the cause of the violent instabilities seen 
in many previous SPMHD calculations with large density 
gradients. To resolve this, we replace the h^g^b} terms in 
these two equations with 

hgb — ~ {hg T hb^ , (21) 

such that 

= — E ^ - ^9 Bi^iWab (kab) , (22) 

at Mo “ Mb ^ ^ 

and for the SPMHD induction equation, 

1 (^) = -4E""^ (Kb) . (23) 

at \ Pa J ha b ^ ^ 

This corrects the instability since in the limit hb hg 
hgb ha and vice versa. The stability seen in Price et al. 
(2012) was simply a product of the larger sink radii pre¬ 
venting the formation of density gradients so extreme that 
this is an issue. Similarly, in the outer regions of the collapse 
simulation where the density gradient is much flatter it does 
not cause a large change in the smoothing length, prevent¬ 
ing an undesirable loss of resolution in these areas. Whilst 
the subtraction in eq. (18) makes the equations of SPMHD 
less conservative; eq. (17) derives entirely from eq. (20) and 
consequently no additional conservation loss is introduced 
by the use of the average h terms. However, the removal of 
the Q terms from eq. (19) and eq. (20) may introduce a very 
small error. 
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Since all other SPH equations (except the den¬ 
sity/smoothing length) contain terms for both ha and hb it 
is not necessary to apply this correction elsewhere. In par¬ 
ticular, we do not apply it to the hyperbolic cleaning terms. 
Even though the density is evaluated using only one smooth¬ 
ing length, this does not contribute to the instability since, 
in effect, p and h are the same parameter and are therefore 
consistent. We did investigate only applying the average h 
method to only one of eq. (22) or eq. (23), and observed that 
this was less stable that applying it to both - which is the 
expected result given that hah will tend towards the larger 
value. 

We also considered a slightly simpler scheme, whereby 
we imposed a minimum smoothing length, Zimin on the whole 
simulation by modifying eq. (10) to be 

ha{ra) = v (~^ + (24) 

which has the desirable properties that it does not intro¬ 
duce any extra loss of conservation and does not change the 
Q terms (since there is no spatial dependence to Zimin). How¬ 
ever, it is difficult to determine a ‘correct value’ of Zimin a 
priori. In addition this formalism would cause benign par¬ 
ticle pairing (see Price (2012) for details) due to particles 
having too many neighbours as well as needlessly sacrific¬ 
ing resolution on other SPH equations. We found that the 
only effective values of Zimin were so large that the loss of 
resolution caused by pairing was in itself a serious problem. 

We also investigated, unsuccessfully, using the average 
of two smoothing kernels, i.e. 

Wab = ^ (Wab {ha) + Wab {hb)) , (25) 

which was little different to the status quo. In the limit where 
ha » hb then Wab(ha) » Wab{hb)- This will result in the 
average in eq. (25) essentially becoming ^Wab{ha)- Whilst 
this would be desirable for one of the two MHD equations, it 
will reduce the smoothing applied to the other substantially. 
Consequently this approach was rejected. 


3 INITIAL CONDITIONS 

The initial conditions for our calculations of protostellar col¬ 
lapse are broadly the same as those in Price et al. (2012). 
However, we use more SPH particles and smaller accre¬ 
tion radii for our sink particles. We begin with a 1.5 mil¬ 
lion SPH particle uniform density sphere of cold gas, more 
than sufficient to resolve a Jeans length according to the 
criteria in Bate & Burkert (1997), placed in a periodic 
box and surrounded by an external medium of ca. 500,000 
warm gas particles. There is a density ratio of 30:1 be¬ 
tween the warm outer medium and the cool sphere with 
a pressure equilibrium between the sphere and the medium. 
Particles are initially laid out on a cubic lattice, the ini¬ 
tial radius of the sphere is rdoud = 4 x 10^^ cm with a 
mass of M = IMq giving an initial density in the sphere 
po = 7.4 X 10~^® g cm“^ The sphere has an initial isother¬ 
mal sound speed Cs = 2.2 x 10^ cm s“^ and we use the 
barotropic equation of state 


P = ci{ 


P 

Pc,l 

Pc,l 


Pc 


,1 



Pc,l < p ^ Pc,2 ^26) 


P > pc,2 


where the two critical densities are given by pc,i = 
10“^^ g cm“^ and pc ,2 = 10“^° g cm“^. This is similar 
to that used, for example, in Machida et al. (2008) with the 
removal of the final 7 = | step at the highest densities. The 
sphere has an initial temperature of approximately 10 K; 
since the outer medium also begins with the same initial 
pressure it has a correspondingly higher initial temperature 
of approximately 300 K. The sphere is set in solid body ro¬ 
tation at Q = 1.77 X 10“^^ rad s“^, such that the magnitude 
of the ratio of rotational to gravitational energy is ~ 0.005, 
within the range observed by Goodman et al. (1993). 


We then define a new parameter, -d, which is the angle 
between the rotation axis of the sphere (which is always 
aligned with the z-axis) and the initial magnetic field. The 
magnetic field is then initially 

Bx = Ho sin P , (27) 

Bz — Ho cos P , (28) 


i.e. when p = 0° the field is aligned with the z-axis. The 
initial magnetic field Ho is determined using the parameter 
p, which is (Nakano & Nakamura 1978; Mac Low & Klessen 
2004) the ratio between the sphere’s mass-flux ratio and the 
critical mass-flux ratio for a spherical cloud, i.e. 


P = 


Z^cloud 

Pcrit 


where. 


(29) 


Z^cloud — 


M 


.So’ 


Pcrit — 


2ci 


TvGpo 


(30) 


with the ratio between the minimum self-collapsing grav¬ 
itational mass obtained from the virial theorem and that 
required for a magnetised astrophysical cloud, ci = 0.53 
as obtained numerically by Mouschovias Sz Spitzer (1976). 
Throughout this paper we use p = 5. 


Sink particles are added to the simulation once a criti¬ 
cal density of 10“^° g cm“^ is achieved, and the usual tests 
are passed (Bate et al. 1995). Our more stable formalism 
allows for smaller sink sizes that used in Price et al. (2012); 
we use an accretion radius of 1 AU as a compromise between 
capturing physics and numerical efficiency. Tests have been 
performed with sink particles with accretion radii of 0.1 AU 
and 0.01 AU. Our sink particle will accrete unconditionally 
once a particle crosses its accretion radius; since all our sim¬ 
ulations are of a collapsing core this should not result in 
any deleterious effects. As in previous work, the sink par¬ 
ticle does not carry a magnetic field - when a particle is 
eliminated from the simulation, the mass is added to the 
sink (which does not exert a hydrodynamic pressure). 
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Normal T = 3.5 T = 4.5 T=5 





-4 -2 0 


log density 


Figure 1. Density cross-sections in the x — y plane for a dif¬ 
ferentially rotating disc. The top row is the ‘standard’ SPMHD 
formalism whilst the bottom is the average h method. The un¬ 
physical bubble caused by the instability discussed above can be 
clearly seen at t = 4.5. 



Figure 2. Total energy and total linear momentum for the ‘stan¬ 
dard’ SPMHD method (solid line) and the average-h method 
(dashed line). The rapid increase in the total energy is correlated 
with the unphysical bubble seen in fig. 1. 


4 ALGORITHM TESTS 

4.1 Isothermal cylinder in a box 

In addition to performing full-scale models, we performed a 
series of low resolution tests. An 8,000 SPH particle isother¬ 
mal cylinder was created inside a periodic box, with a central 
potential provided by a sink particle with a mass 10 times 
that of the cylinder. We use a sink particle rather than a 
potential well to prevent a high density region centred on 
the origin requiring a very short time-step. The inner and 
outer cylinder radii were 0.5 and 5 code units respectively 
and the cylinder was 2.5 units thick (i.e. a height to radius 
ratio of |), and the sink particle had an accretion radius of 
0.3 code units. The cylinder was then set in differential rota- 



Figure 3. Column density plots for the low resolution test with 
the unmodified scheme (top) and our modified average-h scheme 
(bottom) at t = 24940 yrs. The large unphysical explosion can 
be clearly seen in the upper plot, where a large asymmetrical 
bubble of material has been ejected at very high velocities. In 
comparison, the modified scheme forms a collimated jet correctly 
(albeit underresolved due to the low resolution). 


tion, with a velocity profile. The initial velocity was set 
such that the rotation period at 1 unit distance was T = 2. 

A uniform initial sound speed of 0.1 code units was set 
with an isothermal equation of state P{p) = ^up; with an 
initial ratio between hydrodynamic and magnetic pressure 
(the ‘plasma /3’) of /3 ^ 8.4 (corresponding to a nominal 
mass-flux ratio of 5, though this is not a useful measure for a 
differentially rotating cylinder where the effect of magnetic 
braking is much more dominant than magnetic pressure). 
The system was then allowed to evolve. In a correct model, 
we would expect the cylinder to pile up, with material from 
within 1 unit radius moving outwards and more distant ma¬ 
terial spiraling inwards. Some material will fall towards the 
central sink (both due to magnetic and viscous braking) and 
be accreted. The cylinder will also flatten due to rotational 
forces and self gravity, further increasing the density. 

Figure 1 shows cross-sections in the x — y plane for the 
normal SPMHD formalism and our modified one. A clearly 
unphysical bubble like structure can be seen for the nor¬ 
mal code which is not present in our modified method. In 
the original method, the total energy increases rapidly and 
quickly becomes positive; in contrast the average h method 
results in a monotonically decreasing energy, as expected 
for an isothermal equation of state. Earlier, we noted that 
this formalism could, in principle, exhibit poorer momen¬ 
tum conservation than the standard form, however, this is 
in practice undetectable, as seen in fig. 2. 


4.2 Low resolution spherical collapse 

Using the initial conditions discussed more fully in section 3, 
we then performed a low resolution (150,000 SPH particles) 
comparison of the collapse of a spinning magnetised cloud 
core with ^ = 0° using both schemes. Figure 3 shows the 
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Figure 4. Vertical velocity of SPH particles as a function of 
height for both schemes at the same time as in fig. 3. Unlike the 
modified scheme, the original method does not produce a sym¬ 
metrical outfiow at 5 km s“^ (note that this is slower than that 
seen at higher resolutions later since the accretion region is under¬ 
resolved), and the maximum velocities seen are over 40 km s“^. 
Note that the scale on the upper panel is different to that on the 
lower panel. 


situation shortly after the insertion of a sink particle. The 
violent and unphysical explosion in the unmodified code can 
be clearly contrasted with the symmetric bipolar outflow in 
the modified code. 

In both cases, the tests for insertion of a sink particle are 
passed at approximately one free-fall time (tff = 24430 yrs) 
and the explosion happens at between 1.01 tff and 1.02 tff, 
i.e. soon after the critical density for sink creation is reached. 

In the unmodified scheme, a high velocity bubble of ma¬ 
terial is produced and ejected, similar to that seen in fig. 1, 
but in this case the most significant effect is to eject material 
perpendicular to the plane of the disc. This is probably due 
to this being both the rotation and magnetic field axis, and 
therefore the preferred direction for momentum transport 
(similar to how the collimated jet is produced in the modi¬ 
fied scheme). In fig. 4 the velocity in the vertical as a function 
of height is shown, which demonstrates the symmetrical and 
collimated nature of the outflow in the modified code, and 
the much higher and broadly distributed velocities in the 
original code. The collimated jet produced by the average 
h method remains stable until all material has either been 
accreted or ejected and the jet is extinguished. 


5 RESULTS 

5.1 Nature of outflows and jets 

We performed calculations with values of = 

[0°, 10°, 20°, 45°, 60°, 90°]. Figures 5 and 6 show the time 
evolution for these six angles. We note that the the results 


for the 19 = 0° case are broadly the same as in Price et al. 
(2012), albeit with a slightly faster jet velocity - in this case 
~ 8 km s“^. This is expected, since the smaller accretion 
radii used here will allow a faster velocity near the sink 
particle, and since the axial velocity of a collimated jet is 
proportional to the velocity of the matter spiraling in to 
create it this naturally leads to a faster jet (Price et al. 
2003). This is the only significant difference between this 
result and the earlier calculation that used a 5 AU sink, 
showing that our modification to the SPH equations does 
not cause numerical artifacts or errors of its own. 

The most striking result is the lack of any real outflow 
at all from the i9 = [60,90]° models. Whilst all shallower 
angles produce an outflow of some significance, this only 
takes the form of a collimated jet for i9 ^ 20°, and can 
only be sustained when i9 ^ 10°. Similarly, the pseudo-disc 
which is clearly defined for i9 = 0° is either disrupted or, in 
the most misaligned cases, does not form at all. Since the 
formation of a stable bipolar outflow requires a stable and 
defined disc structure, this naturally prevents a substantial 
outflow being formed. In the intermediate i9 = 45° case, the 
pseudo-disc formed is highly disrupted but still manages to 
drive a broad, albeit slower, outflow. 

We observe in fig. 8 that as the magnetic field geometry 
near the sink becomes very complicated bubbles of material 
driven by magnetic pressure form and disrupt the accretion 
of matter into a disc. For i9 ^ 60° this is sufficient to sup¬ 
press the formation of a disc and outflow altogether, whilst 
for shallower angles the outflow simply becomes more diffuse 
and less structured. Compared to the aligned case, even set- 
ing i9 ^ 10° has an effect on the pseudo-disc (hg. 9). In and 
of itself, such a structure should not prevent a collimated jet 
being produced - and indeed one is seen in both the 10 and 
20° models with a similar velocity to the simpler aligned 
case. Figure 7 shows the evolution of the jet for i9 = 10°. 
Initially, the most central region of the pseudo-disc aligns 
perpendicularly with the rotation axis, and consequently a 
jet is produced parallel to this axis, whilst the outer regions 
remain aligned perpendicular to the held axis. As the system 
evolves, the central portion warps, this causes the outflow 
to realign with the held axis, causing the kink seen in hg. 5 
at the extremities of the jet. 

In hg. 10 we plot a series of density cross-sections for 
the i9 = 20° model. The same disc warping seen at 10° is 
present, however, unlike the shallower angle the jet produced 
here is disrupted within 700 years of forming. The outflow 
continues, as seen in hg. 5, and has a prohle broadly similar 
to shallower angles but without a central jet. This loss of 
collimation appears to be the result of the formation of a 
bubble of material near the protostar which pushes material 
away from the core and thereby prevent the formation of a 
high velocity region of the pseudo-disc which can collimate 
an outflow. This is not a true accretion disc with a Keple- 
rian velocity prohle and is simply a result of rotational and 
magnetic forces forcing material into a disc-like structure. 
Whilst magnetic bubbles of this nature have been seen be¬ 
fore (see, e.g. Zhao et al. (2011), Krasnopolsky et al. (2012)) 
it is unclear whether this is a real effect or due to the lack 
of any physical resistivity or other non-ideal MHD effects 
(e.g. ambipolar diffusion or the Hall effect) in our model. 
Resistivity may allow for magnetic reconnection and hence 
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Figure 5. Column density projection plots for -d = [0°, 10°, 20°] (across the page) at 4 different times (down the page). At these shallow 
angles, a prominent collimated outflow aligned with the field axis is always formed, however in the 7? = 20° case this is eventually 
disrupted and becomes puffy and uncollimated hy t = 27870 years. The rotation axis is along the z-axis. 
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Figure 6. Column density projection plots for 'd = [45°, 60°, 90°] (across the page) at 4 different times (down the page). Note that the 
scale here is different to fig. 5 to show the more complicated inner structures. At these substantially steeper angles no collimated outflow 
is produced at all and for -d > 45° the outflow is very heavily suppressed. 
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Figure 7. Cross-sections of density in the z — x plane for -d = 10°. 
Whilst the outer regions of the pseudo-disc align perpendicular 
to the magnetic field, the innermost region exhibits a more com¬ 
plicated structure. As it deforms, the collimated jet changes from 
being parallel to the rotation axis to being parallel to the field 
axis. 



Figure 8. Cross-sections of |J| at t = 25420 yrs for 4 different 
values of t?. The magnetic field geometry is significantly more 
complicated in the latter two cases, and this corresponds with a 
substantially changed outfiow. 

the conversion of magnetic energy into thermal energy and 
may act to stabilise the accretion structures in this case. 

A similar effect can be seen for 'd = 45°, however in 
this case the effect is so pronounced that the outflow is com¬ 
pletely disrupted and very large, broadly symmetrical ‘puffy’ 
outflow can be seen. Comparing the shallower angles with 
this model in fig. 11 shows that the actual velocity profile 
consists of two zones: a narrow region just above and below 
the sink that is aligned with the rotation axis, and a diffuse 
yet still fast zone further out which broadly follows the field 
axis. The alignment of the outflow initially with the rota- 



-18 -16 -14 

log density [g/cnn^] 


Figure 9. Cross-sections of density in the x-y plane (left column) 
and z-y plane (right column) at t = 27870 yrs for the same values 
of 'd as in fig. 8. The uniform pseudo-disc, perpendicular to the 
rotation axis, can be clearly seen in the aligned case, as can the 
formation of a bar like structure for 'd = 10°. The jets in the lower 
two pairs of plots are, as expected given the complicated velocity 
profile, completely disrupted. 

tion axis is caused by the central region of the pseudo-disc 
being aligned (as shown in fig. 9) perpendicular to the rota¬ 
tion axis. As the outflow moves away, it will be acted on by 
magnetic force which will re-orient it to align with the field 
axis so that the fluid moves along the field lines. Unlike for 
'd — 20° and 60°, the bubble structures seen at this angle 
remain generally symmetrical; a likely explanation for this 
is simply that for 45° the magnitude of the field is identi¬ 
cal in both the x and z direction and consequently there is 
nothing to create an asymmetry. 


© 2015 RAS, MNRAS 000, 1-13 
















10 B. T. Lewis, M. R. Bate, and D. J. Price 



Figure 10. Density cross-sections for 'd = 20°. The effect of the 
complicated field geometry in expelling material from around the 
protostar, and thus distrupting the collimated jet, can be clearly 
seen. 


We note that the jets seen in the more aligned cases 
(ca. 8 km s“^ ) are significantly slower than those observed 
in observations of Herbig-Haro objects (where velocities of 
>100 km s~^ are commonly seen) and are also slower than 
the fastest velocities observed in Bate et al. (2014). This is 
caused by the relatively large sink radius which is orders of 
magnitude larger than the stellar core, and the consequent 
relatively slow maximum Keplerian velocity (in this case at 
1 AU). Since jet velocity is proportional to the velocity of the 
antecedent accretion disc (Price et al. 2003) this naturally 
produces a slower jet. 

The lack of any outflow, as opposed to a simple puffing 
out of the core region, in the ^ = 90° model is expected since 
a field geometry this extreme will naturally lead to a very 
complicated distribution of field and mass near the sink. In 
fig. 12 we compare the plasma /3 in the aligned and most 
non-aligned cases. This clearly shows that bubbles of mate¬ 
rial, driven by magnetic pressure disrupt the accretion disc 
whilst in the aligned case this hydrodynamically dominated 
disc persists all the way to the centre and ultimately the 
accretion region. As seen previously, fig. 11 shows that the 
velocities in this case are substantially lower than for the 
aligned model (the maximum outflow velocity is approxi¬ 
mately 2 km s“^, compared to 8 km s“^ when = 0°) pro¬ 
viding further evidence that magnetic pressure is influencing 
the collapse. 

We obtain generally similar morphologies to those seen 
in Ciardi & Hennebelle (2010), where the outflows become 
increasingly disrupted as 'd is increased. For example, at= 
45° a puffy outflow with no distinct jet is obtained in both 
models, and at 27000 yrs this outflow is 2000 — 3000 
AU in size; and at higher angles the outflow is suppressed 
until being essentially extinguished at = 90°. Both models 
also agree for = 20° until ~ 22000 yrs when Ciardi & 
Hennebelle (2010) stop, however, we observe that the jet is 
subsequently disrupted and replaced with a diffuse outflow. 


Figure 11. Cross-sections of |v| at T = 27870 yrs (i.e. late enough 
that the magnetic effects seen in fig. 8 will have been able to dis¬ 
rupt the outflow). The two upper plots show the strongly colli¬ 
mated outflows seen for low values of 7? whilst the 7 ? = 45 and 
90° plots show the more complicated, diffuse outflows which at 
steeper angles essentially cease altogether. 



Figure 12. Plasma /3 cross-sections for 7 ? = 0° and 90°. In the 
aligned case, a pressure supported pseudo-disc (aligned perpen¬ 
dicularly to the field) is present, whilst for the misaligned case 
it is completely disrupted by magnetic pressure several hundred 
AU away. 

5.2 Accretion 

In the early part of the simulation, after the sink is inserted, 
we observe very rapid accretion. This rate decreases over 
time, both due to matter being expelled from the core by 
outflows and also due to the dynamics of the collapse (Foster 
& Chevalier 1993), however, the eventual accretion rate does 
depend on the value of 'd used. 

Figure 13 shows the mass accreted by the sink particle 
for each model as a function of time. We would assume that 
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25000 26000 27000 28000 

Time [yr] 

Figure 13. Accretion of mass by the sink particle for -i? = 0° 
(solid), 10° (dotted), 20° (short dashes), 45° (long dashes), 60° 
(dot - short dashes), and 90° (dot - long dashes). The 10° model 
accretes slightly faster than the aligned model; all steeper angles 
except 90° accrete less material, although 20° initially follows 10°. 
There is a sharp ‘knee’ in the 90° model at approximately 28000 
yrs. 

for ever increasing values of the accretion rate will fall as 
the complicated field effects (ranging from loops of material 
to the more extreme 90° model) push material away from 
the sink particle. Between ~ 20° and ~ 60° the accretion 
rate does indeed fall for steeper angles. However this trend 
markedly reverses for the steepest angles, in particular the 
rate for 90° is substantially faster and is observed to increase 
sharply at a ‘knee’ rather than decrease. This is a counter¬ 
intuitive result given the structure of the cloud core and the 
very low plasma [3 in this regime. This differs substantially 
from the result Ciardi & Hennebelle (2010) who found that 
there was no regime in which the accretion rate fell. Ciardi & 
Hennebelle (2010) developed an analytical model (from the 
model for an isothermal collapsing sphere in Hunter (1977)) 
whereby 

Mcore (t) = Tae Mnf ^1 - CXp , (31) 

where Minf is a constant determined from the sound speed 
of the medium and 

Tae OC ^ . (32) 

COS ly 

Equations (31) and (32) will produce ever faster accretion 
rates as increases without the decrease we see for values 
of 20° ^ < 90° - and clearly breaks down when = 90°. 

This is simply a result of the assumptions made, i.e. that 
the sphere collapses held up only by magnetic pressure (and 
by material being removed by the outflow). We do obtain 
an increased accretion rate in the 10° model (and initially, 
before magnetic bubbles distrupt the pseudo-disc, in the 20° 
model) which lends some support to this hypothesis. In con- 



Figure 14. Latitude-longitude maps of particle accretion for the 
first 4,000 after a sink is inserted. The upper plot is for the = 0° 
case and the lower plot for -d = 90°. Darker colours, which denote 
a greater number of particles accreted show that for the aligned 
case the accretion is clearly disc-linked, and conversely that the 
misaligned accretion covers a substantially larger solid angle. The 
color scale is normalised so that white represents a pixel with no 
accreted particles and black represents the pixel with the highest 
number of accreted particles. 

trast, the complicated bubble structures seen at larger angles 
disrupt the accretion process and therefore cause a reduc¬ 
tion in the core mass. At the largest values of -d, the outflow 
is so suppressed, however, that the accretion rate actually 
increases. 

We see in fig. 14 one of the causes of this increase in 
accretion rate at -d = 90°. For strongly aligned fields and ro¬ 
tation axes, the accretion process can only happen along the 
edge of the disc. This remains true even as the disc itself is 
disrupted by magnetic effects - in essence, rather than being 
a constant equatorial line, material is accreted at the edges 
of the bubbles and other disturbances in the pseudo-disc. For 
much larger angles, the solid angle over which accretion can 
occur is much larger - whilst the rotational forces are trying 
to hold the material into a disc-like structure, the magnetic 
bubbles formed completely disrupt this. There is a still a 
general preference to accrete material along the equator of 
the sink, rather than the poles because the material is still 
spinning. 


6 CONCLUSION 

We have devised a small modification to the equations 
of SPMHD, implementing the induction equation and 
anisotropic force equation using an average smoothing 
length formalism, to eliminate an instability that has hith¬ 
erto prevented work using smaller sink particles and mis¬ 
aligned magnetic fields in SPMHD models of collapsing mag¬ 
netised cloud cores. We produced six models of these cloud 
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cores, with angles between the field and rotation axes rang¬ 
ing from 0° (aligned) to 90°. In each simulation, we were able 
to follow the collapse from an initially cold sphere for over 
26,000 years to the formation of a dense core, accretion disc 
and outflows. 

We observe that the nature of the outflows observed 
varies strongly with the initial field geometry. For the 
aligned, i.e. -d = 0°, model we observe an outflow similar 
to that obtained in previous work (Price et al. 2012) albeit 
with a slightly faster collimated jet due to the smaller sink 
radius, with the characteristic jet and envelope seen in pre¬ 
vious work. For misaligned fields, we are able to divide the 
nature of the outflow into three regimes: for very shallow 
angles (^ 10°) a collimated jet is produced; for moderate 
angles (20° — 45°) an initial collimated jet is rapidly dis¬ 
rupted by the increasingly wound up magnetic field which 
produces large magnetic bubbles and the outflow becomes 
much more diffuse and unstructured. Finally, for the steepest 
angles the outflow is substantially suppressed and becomes 
more spherical. 
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